###### *Immigration and public support for political systems in Europe*
###### Replication code for TSCS analyses
###### Christopher Claassen 2022

library(texreg)
library(dplyr)
library(tidyr)
library(ggplot2)
library(RColorBrewer)
library(viridis)
library(psych)
library(plm)
library(MASS)
library(systemfit)
library(lavaan)

# set up some options
plotcol7 = viridis(7, alpha=1, begin=0.1, end=0.9)

# working directory
WD = rstudioapi::getActiveDocumentContext()$path 
setwd(dirname(WD))
print( getwd() )

# load data
migrat = read.csv("EU_migrat_data_PoP.csv")

# check data
migrat$Country = as.factor(migrat$Country)
levels(migrat$Country)
table(migrat$Year)
summary(migrat)

# create panel dataset
mig.plm = pdata.frame(migrat, index = c("Country", "Year"))
sort(rowSums(table(index(mig.plm), useNA = "ifany")))
pdim(mig.plm) # n=30, T=41, N=1230

# create western europe dataset
cnts.drop = c("Croatia", "Czechia", "Hungary", "Poland", "Romania", "Slovakia", "Slovenia", 
              "Latvia", "Estonia", "Lithuania")
migrat.weur = migrat[!migrat$Country %in% cnts.drop,]
migrat.weur$Country = as.factor(as.character(migrat.weur$Country))
mig.plm.we = pdata.frame(migrat.weur, index = c("Country", "Year"))
sort(rowSums(table(index(mig.plm.we), useNA = "ifany")))
pdim(mig.plm.we) # n=20, T=41, N=820


### Create some lagged and within measures for plotting

# Create lagged measures

cnts = levels(as.factor(as.character(migrat$Country)))
yrs = 1990:2020

for(i in 1:length(cnts)) {
  for(j in yrs) {
    migrat[migrat$Country==cnts[i] & migrat$Year==j, "Immig_m1"] =
      migrat[migrat$Country==cnts[i] & migrat$Year==j-1, "Immig_pc"]
    migrat[migrat$Country==cnts[i] & migrat$Year==j, "Forborn_m1"] =
      migrat[migrat$Country==cnts[i] & migrat$Year==j-1, "Forborn_stck_pc"]
  }
}

# create within and between measures

mig.trim = migrat
mig.trim = mig.trim[complete.cases(mig.trim[, c("Trust", "Forborn_m1", "Immig_m1")]), ]

mig.trim$Forborn_stck_btw = ave(mig.trim$Forborn_stck_pc, mig.trim$Country, FUN=function(x) mean(x, na.rm=TRUE))
mig.trim$Forborn_stck_within = mig.trim$Forborn_stck_pc - mig.trim$Forborn_stck_btw
mig.trim$Immig_btw = ave(mig.trim$Immig_pc, mig.trim$Country, FUN=function(x) mean(x, na.rm=TRUE))
mig.trim$Immig_within = mig.trim$Immig_pc - mig.trim$Immig_btw


### Time series plots


cnts = levels(migrat$Country)

# trust and immigration

pdf("Fig_1-trust_immig_ts.pdf", height=7.5, width=6.5)
par(mfrow=c(6,5), mar=c(1.2,1,0.2,0.2), tcl=-0.15, cex=0.9)
for(i in 1:length(cnts)) {
  plot(x=migrat[migrat$Country==cnts[i], "Year"],
       y=migrat[migrat$Country==cnts[i], "Trust"],
       type="n", xlim=c(1980, 2020), ylim=c(-3, 4.2), xlab="", ylab="", main="", axes=FALSE)
  grid(col=rgb(0,0,0,0.5), lwd=0.5)
  lines(x=migrat[migrat$Country==cnts[i], "Year"],
        y=migrat[migrat$Country==cnts[i], "Immig_pc"],
        col=plotcol7[6], lwd=2)
  lines(x=migrat[migrat$Country==cnts[i], "Year"],
        y=migrat[migrat$Country==cnts[i], "Trust"],
        col=plotcol7[3], lwd=2)
  axis(side=1, at=c(1980, 1990, 2000, 2010, 2020), labels=c(1980, 1990, 2000, 2010, 2020), mgp=c(1, -0.2, 0), 
       cex.axis=0.65)
  axis(side=2, at=c(-2,0,2,4), mgp=c(1,0.1,0), cex.axis=0.65)
  text(cnts[i], x=1980, y=3.8, cex=0.75, adj=0)
  box()
}
text(x=2023, y=1.3, labels="Immigration rate (% of pop)", pos=2, col=plotcol7[6], cex=0.6)
text(x=2023, y=-0.9, labels="Institutional trust", pos=2, col=plotcol7[3], cex=0.6)
dev.off()


# satisfaction and immigration

pdf("Fig_S5-satis_immig_ts.pdf", height=7.5, width=6.5)
par(mfrow=c(6,5), mar=c(1.2,1,0.2,0.2), tcl=-0.15, cex=0.9)
for(i in 1:length(cnts)) {
  plot(x=migrat[migrat$Country==cnts[i], "Year"],
       y=migrat[migrat$Country==cnts[i], "Satis"],
       type="n", xlim=c(1980, 2020), ylim=c(-3, 4.2), xlab="", ylab="", main="", axes=FALSE)
  grid(col=rgb(0,0,0,0.5), lwd=0.5)
  lines(x=migrat[migrat$Country==cnts[i], "Year"],
        y=migrat[migrat$Country==cnts[i], "Immig_pc"],
        col=plotcol7[6], lwd=2)
  lines(x=migrat[migrat$Country==cnts[i], "Year"],
        y=migrat[migrat$Country==cnts[i], "Satis"],
        col=plotcol7[2], lwd=2)
  axis(side=1, at=c(1980, 1990, 2000, 2010, 2020), labels=c(1980, 1990, 2000, 2010, 2020), mgp=c(1, -0.2, 0), 
       cex.axis=0.65)
  axis(side=2, at=c(-2,0,2,4), mgp=c(1,0.1,0), cex.axis=0.65)
  text(cnts[i], x=1980, y=3.8, cex=0.75, adj=0)
  box()
}
text(x=2023, y=1, labels="Immigration rate (% of pop)", pos=2, col=plotcol7[6], cex=0.6)
text(x=2024, y=-1, labels="Satisfaction with democracy", pos=2, col=plotcol7[2], cex=0.6)
dev.off()


# democratic support and immigration

pdf("Fig_S6-mood_immig_ts.pdf", height=7.5, width=6.5)
par(mfrow=c(6,5), mar=c(1.2,1,0.2,0.2), tcl=-0.15, cex=0.9)
for(i in 1:length(cnts)) {
  plot(x=migrat[migrat$Country==cnts[i], "Year"],
       y=migrat[migrat$Country==cnts[i], "SupDem"],
       type="n", xlim=c(1980, 2020), ylim=c(-3, 4.2), xlab="", ylab="", main="", axes=FALSE)
  grid(col=rgb(0,0,0,0.5), lwd=0.5)
  lines(x=migrat[migrat$Country==cnts[i], "Year"],
        y=migrat[migrat$Country==cnts[i], "Immig_pc"],
        col=plotcol7[6], lwd=2)
  lines(x=migrat[migrat$Country==cnts[i], "Year"],
        y=migrat[migrat$Country==cnts[i], "SupDem"],
        col=plotcol7[4], lwd=2)
  axis(side=1, at=c(1980, 1990, 2000, 2010, 2020), labels=c(1980, 1990, 2000, 2010, 2020), mgp=c(1, -0.2, 0), 
       cex.axis=0.65)
  axis(side=2, at=c(-2,0,2,4), mgp=c(1,0.1,0), cex.axis=0.65)
  text(cnts[i], x=1980, y=3.8, cex=0.75, adj=0)
  box()
}
text(x=2023, y=1, labels="Immigration rate (% of pop)", pos=2, col=plotcol7[6], cex=0.6)
text(x=2024, y=-1, labels="Support for democracy", pos=2, col=plotcol7[4], cex=0.6)
dev.off()



### TSCS Models

## Trust as DV, dynamic FE

form1 = diff(Trust, lag=1) ~ plm::lag(Trust, 1:2) + plm::lag(Forborn_stck_pc, 1) + plm::lag(Immig_pc, 1) 
form2 = diff(Trust, lag=1) ~ plm::lag(Trust, 1:2) + diff(GDP_grwth_WB, lag=1) + plm::lag(GDP_grwth_WB, 1) + 
  plm::lag(Unempl_WB, 1) + plm::lag(Corrup_VD, 1) + plm::lag(Libdem_VD, 1) + plm::lag(Forborn_stck_pc, 1) + 
  plm::lag(Immig_pc, 1)
form3 = update.formula(form2, . ~ . + plm::lag(FRight_seat_pr, 1) + plm::lag(BMI_i_ext, 1) )
form4 = update.formula(form2, . ~ . - plm::lag(Immig_pc, 1) + plm::lag(Imm_mus_pc, 1) )
form5 = update.formula(form2, . ~ . - plm::lag(Immig_pc, 1) + plm::lag(Imm_nonEU_pc, 1) )
form6 = update.formula(form2, . ~ . + plm::lag(Imm_sal_2, 0) )

trust.mods.ip = lapply(list(form1, form2, form3, form4, form5, form6), 
                       plm, data=mig.plm, model="within", effect="indiv")
sapply(trust.mods.ip, function(x) round(pwartest(x, order=1)[["p.value"]][[1]], 3) )
sapply(trust.mods.ip, function(x) round(pcdtest(x, test="cd")[["p.value"]][[1]], 3) )
sapply(trust.mods.ip, function(x) round(sigma(x), 3) )

vcm = lapply(trust.mods.ip, function(x) plm::vcovSCC(x, cluster="time")) # cross-sectional heteroskedasticity and correlation
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(trust.mods.ip)) {pval[[i]] = pt(abs(coef(trust.mods.ip[[i]]) / rse[[i]]), df=400, lower.tail=FALSE) * 2}
screenreg(trust.mods.ip, override.se=rse, override.pvalues=pval, digits=3, single.row=FALSE, leading.zero=FALSE, 
          stars=c(0.05))

vars = c("Institutional trust t-1", "Institutional trust t-2", "% foreign-born t-1", "Immigration rate t-1", 
         "Delta GDP growth per capita t0", "GDP growth per capita t-1", "Unemployment rate t-1", 
         "Corruption t-1", "Liberal democracy t-1", "Far right seat share t-1", "Immigrant integration policy t-1", 
         "Muslim immigration rate t-1", "Non-EU immigration rate", "Concern about immigration t0")
htmlreg(trust.mods.ip, file="table1.htm", override.se=rse, override.pvalues=pval, digits=3,
        single.row=TRUE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
        caption="Immigration and institutional trust, 1-way FE mods, DKSEs")
texreg(trust.mods.ip, file="table1.txt", override.se=rse, override.pvalues=pval, digits=3,
       single.row=FALSE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
       caption="Immigration and institutional trust, 1-way FE mods, DKSEs",
       booktabs=TRUE, dcolumn=TRUE,)


## Satis as DV, dynamic FE

form1 = diff(Satis, lag=1) ~ plm::lag(Satis, 1:2) + plm::lag(Forborn_stck_pc, 1) + plm::lag(Immig_pc, 1) 
form2 = diff(Satis, lag=1) ~ plm::lag(Satis, 1:2) + diff(GDP_grwth_WB, lag=1) + plm::lag(GDP_grwth_WB, 1) + 
  plm::lag(Unempl_WB, 1) + plm::lag(Corrup_VD, 1) + plm::lag(Libdem_VD, 1) + plm::lag(Forborn_stck_pc, 1) + 
  plm::lag(Immig_pc, 1)
form3 = update.formula(form2, . ~ . + plm::lag(FRight_seat_pr, 1) + plm::lag(BMI_i_ext, 1) )
form4 = update.formula(form2, . ~ . - plm::lag(Immig_pc, 1) + plm::lag(Imm_mus_pc, 1) )
form5 = update.formula(form2, . ~ . - plm::lag(Immig_pc, 1) + plm::lag(Imm_nonEU_pc, 1) )
form6 = update.formula(form2, . ~ . + plm::lag(Imm_sal_2, 0) )

satis.mods.ip = lapply(list(form1, form2, form3, form4, form5, form6), 
                       plm, data=mig.plm, model="within", effect="indiv")
sapply(satis.mods.ip, function(x) round(pwartest(x, order=1)[["p.value"]][[1]], 3) )
sapply(satis.mods.ip, function(x) round(pcdtest(x, test="cd")[["p.value"]][[1]], 3) )
sapply(satis.mods.ip, function(x) round(sigma(x), 3) )
vcm = lapply(satis.mods.ip, function(x) plm::vcovSCC(x, cluster="time")) # cross-sectional heteroskedasticity and correlation
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(satis.mods.ip)) {pval[[i]] = pt(abs(coef(satis.mods.ip[[i]]) / rse[[i]]), df=400, lower.tail=FALSE) * 2}
screenreg(satis.mods.ip, override.se=rse, override.pvalues=pval, digits=3, single.row=FALSE, leading.zero=FALSE, 
          stars=c(0.05))

vars = c("Satisfaction with democracy t-1", "Satisfaction with democracy t-2", "% foreign-born t-1", "Immigration rate t-1", 
         "Delta GDP growth per capita t0", "GDP growth per capita t-1", "Unemployment rate t-1", 
         "Corruption t-1", "Liberal democracy t-1", "Far right seat share t-1", "Immigrant integration policy t-1", 
         "Muslim immigration rate t-1", "Non-EU immigration rate", "Concern about immigration t0")
htmlreg(satis.mods.ip, file="table2.htm", override.se=rse, override.pvalues=pval, digits=3, 
        single.row=TRUE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
        caption="Immigration and satisfaction with democracy, 1-way FE mods, DKSEs")
texreg(satis.mods.ip, file="table2.txt", override.se=rse, override.pvalues=pval, digits=3,
       single.row=FALSE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
       caption="Immigration and satisfaction with democracy, 1-way FE mods, DKSEs",
       booktabs=TRUE, dcolumn=TRUE,)


## Democratic mood as DV, dynamic FE

form1 = diff(SupDem, lag=1) ~ plm::lag(SupDem, 1:2) + plm::lag(Forborn_stck_pc, 1) + plm::lag(Immig_pc, 1) 
form2 = diff(SupDem, lag=1) ~ plm::lag(SupDem, 1:2) + diff(GDP_grwth_WB, lag=1) + plm::lag(GDP_grwth_WB, 1) + 
  plm::lag(Unempl_WB, 1) + plm::lag(Corrup_VD, 1) + plm::lag(Libdem_VD, 1) + plm::lag(Forborn_stck_pc, 1) + 
  plm::lag(Immig_pc, 1)
form3 = update.formula(form2, . ~ . + plm::lag(FRight_seat_pr, 1) + plm::lag(BMI_i_ext, 1) )
form4 = update.formula(form2, . ~ . - plm::lag(Immig_pc, 1) + plm::lag(Imm_mus_pc, 1) )
form5 = update.formula(form2, . ~ . - plm::lag(Immig_pc, 1) + plm::lag(Imm_nonEU_pc, 1) )
form6 = update.formula(form2, . ~ . + plm::lag(Imm_sal_2, 0) )

mood.mods.ip = lapply(list(form1, form2, form3, form4, form5, form6), 
                       plm, data=mig.plm, model="within", effect="indiv")
sapply(mood.mods.ip, function(x) round(pwartest(x, order=1)[["p.value"]][[1]], 3) )
sapply(mood.mods.ip, function(x) round(pcdtest(x, test="cd")[["p.value"]][[1]], 3) )
sapply(mood.mods.ip, function(x) round(sigma(x), 3) )
vcm = lapply(mood.mods.ip, function(x) plm::vcovSCC(x, cluster="time")) # cross-sectional heteroskedasticity and correlation
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(mood.mods.ip)) {pval[[i]] = pt(abs(coef(mood.mods.ip[[i]]) / rse[[i]]), df=400, lower.tail=FALSE) * 2}
screenreg(mood.mods.ip, override.se=rse, override.pvalues=pval, digits=3, single.row=FALSE, leading.zero=FALSE, 
          stars=c(0.05))

vars = c("Democratic support t-1", "Democratic support t-2", "% foreign-born t-1", "Immigration rate t-1", 
         "Delta GDP growth per capita t0", "GDP growth per capita t-1", "Unemployment rate t-1", 
         "Corruption t-1", "Liberal democracy t-1", "Far right seat share t-1", "Immigrant integration policy t-1", 
         "Muslim immigration rate t-1", "Non-EU immigration rate", "Concern about immigration t0")
htmlreg(mood.mods.ip, file="table3.htm", override.se=rse, override.pvalues=pval, digits=3, 
        single.row=TRUE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
        caption="Immigration and support for democracy, 1-way FE mods, DKSEs")
texreg(mood.mods.ip, file="table3.txt", override.se=rse, override.pvalues=pval, digits=3,
       single.row=FALSE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
       caption="Immigration and democratic support, 1-way FE mods, DKSEs",
       booktabs=TRUE, dcolumn=TRUE,)



### ESS Models

## ESS trust

form0 = ess_trust_all ~ diff(GDP_grwth_WB, lag=1) + plm::lag(GDP_grwth_WB, 1) + plm::lag(Unempl_WB, 1) + 
  plm::lag(Corrup_VD, 1) + plm::lag(Libdem_VD, 1) 
form1 = update.formula(form0, . ~ . + plm::lag(Forborn_stck_pc, 1) + plm::lag(Immig_pc, 1))
form2 = update.formula(form0, . ~ . + plm::lag(Forborn_stck_pc, 1) + plm::lag(Imm_mus_pc, 1) )
form3 = update.formula(form1, . ~ . + plm::lag(ess_imm_sup_all, 0) + plm::lag(ess_satlife_all, 0)  )

form0 = ess_trust_nbrn ~ diff(GDP_grwth_WB, lag=1) + plm::lag(GDP_grwth_WB, 1) + plm::lag(Unempl_WB, 1) + 
  plm::lag(Corrup_VD, 1) + plm::lag(Libdem_VD, 1) 
form4 = update.formula(form0, . ~ . + plm::lag(Forborn_stck_pc, 1) + plm::lag(Immig_pc, 1))
form5 = update.formula(form0, . ~ . + plm::lag(Forborn_stck_pc, 1) + plm::lag(Imm_mus_pc, 1) )
form6 = update.formula(form4, . ~ . + plm::lag(ess_imm_sup_nbrn, 0) + plm::lag(ess_satlife_nbrn, 0)  )

trust.mods.ess = lapply(list(form1, form2, form3, form4, form5, form6), 
                        plm, data=mig.plm, model="within", effect="indiv")
sapply(trust.mods.ess, function(x) round(pwartest(x, order=1)[["p.value"]][[1]], 3) )
sapply(trust.mods.ess, function(x) round(pcdtest(x, test="cd")[["p.value"]][[1]], 3) )
sapply(trust.mods.ess, function(x) round(sigma(x), 3) )
vcm = lapply(trust.mods.ess, function(x) plm::vcovBK(x, cluster="group")) # heteroskedasticity and serial correlation
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(trust.mods.ess)) {pval[[i]] = pt(abs(coef(trust.mods.ess[[i]]) / rse[[i]]), df=200, lower.tail=FALSE) * 2}
screenreg(trust.mods.ess, override.se=rse, override.pvalues=pval, digits=3, single.row=FALSE, leading.zero=FALSE, 
          stars=c(0.05))

vars = c("Delta GDP growth per capita t0", "GDP growth per capita t-1", "Unemployment rate t-1", "Corruption t-1", 
         "Liberal democracy t-1", "% foreign-born t-1", "Immigration rate t-1", "Muslim immigration rate t-1",
         "Evaluations of immigration t0", "Life satisfaction t0", "Evaluations of immigration t0", "Life satisfaction t0")
htmlreg(trust.mods.ess, file="table4.htm", override.se=rse, override.pvalues=pval, digits=3, 
        single.row=TRUE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
        caption="Immigration and institutional trust, ESS only, 1-way FE mods, DKSEs")
texreg(trust.mods.ess, file="table4.txt", override.se=rse, override.pvalues=pval, digits=3,
       single.row=FALSE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
       caption="Immigration and institutional trust, ESS only, 1-way FE mods, DKSEs",
       booktabs=TRUE, dcolumn=TRUE,)


## ESS satisfaction

form0 = ess_satis_all ~ diff(GDP_grwth_WB, lag=1) + plm::lag(GDP_grwth_WB, 1) + plm::lag(Unempl_WB, 1) + 
  plm::lag(Corrup_VD, 1) + plm::lag(Libdem_VD, 1) 
form1 = update.formula(form0, . ~ . + plm::lag(Forborn_stck_pc, 1) + plm::lag(Immig_pc, 1))
form2 = update.formula(form0, . ~ . + plm::lag(Forborn_stck_pc, 1) + plm::lag(Imm_mus_pc, 1) )
form3 = update.formula(form1, . ~ . + plm::lag(ess_imm_sup_all, 0) + plm::lag(ess_satlife_all, 0)  )

form0 = ess_satis_nbrn ~ diff(GDP_grwth_WB, lag=1) + plm::lag(GDP_grwth_WB, 1) + plm::lag(Unempl_WB, 1) + 
  plm::lag(Corrup_VD, 1) + plm::lag(Libdem_VD, 1) 
form4 = update.formula(form0, . ~ . + plm::lag(Forborn_stck_pc, 1) + plm::lag(Immig_pc, 1))
form5 = update.formula(form0, . ~ . + plm::lag(Forborn_stck_pc, 1) + plm::lag(Imm_mus_pc, 1) )
form6 = update.formula(form4, . ~ . + plm::lag(ess_imm_sup_nbrn, 0) + plm::lag(ess_satlife_nbrn, 0)  )

satis.mods.ess = lapply(list(form1, form2, form3, form4, form5, form6), 
                        plm, data=mig.plm, model="within", effect="indiv")
sapply(satis.mods.ess, function(x) round(pwartest(x, order=1)[["p.value"]][[1]], 3) )
sapply(satis.mods.ess, function(x) round(pcdtest(x, test="cd")[["p.value"]][[1]], 3) )
sapply(satis.mods.ess, function(x) round(sigma(x), 3) )
vcm = lapply(satis.mods.ess, function(x) plm::vcovBK(x, cluster="group")) # heteroskedasticity and serial correlation
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(satis.mods.ess)) {pval[[i]] = pt(abs(coef(satis.mods.ess[[i]]) / rse[[i]]), df=200, lower.tail=FALSE) * 2}
screenreg(satis.mods.ess, override.se=rse, override.pvalues=pval, digits=3, single.row=FALSE, leading.zero=FALSE, 
          stars=c(0.05))

vars = c("Delta GDP growth per capita t0", "GDP growth per capita t-1", "Unemployment rate t-1", "Corruption t-1", 
         "Liberal democracy t-1", "% foreign-born t-1", "Immigration rate t-1", "Muslim immigration rate t-1",
         "Evaluations of immigration t0", "Life satisfaction t0", "Evaluations of immigration t0", "Life satisfaction t0")
htmlreg(satis.mods.ess, file="table5.htm", override.se=rse, override.pvalues=pval, digits=3, 
        single.row=TRUE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
        caption="Immigration and satisfaction with democracy, ESS only, 1-way FE mods, DKSEs")
texreg(satis.mods.ess, file="table5.txt", override.se=rse, override.pvalues=pval, digits=3,
       single.row=FALSE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
       caption="Immigration, and satisfaction with democracy, ESS only",
       booktabs=TRUE, dcolumn=TRUE,)


### Western Europe models

## Trust as DV, dynamic FE

form1 = diff(Trust, lag=1) ~ plm::lag(Trust, 1:2) + plm::lag(Forborn_stck_pc, 1) + plm::lag(Immig_pc, 1) 
form2 = diff(Trust, lag=1) ~ plm::lag(Trust, 1:2) + diff(GDP_grwth_WB, lag=1) + plm::lag(GDP_grwth_WB, 1) + 
  plm::lag(Unempl_WB, 1) + plm::lag(Corrup_VD, 1) + plm::lag(Libdem_VD, 1) + plm::lag(Forborn_stck_pc, 1) + 
  plm::lag(Immig_pc, 1)
form3 = update.formula(form2, . ~ . + plm::lag(FRight_seat_pr, 1) + plm::lag(BMI_i_ext, 1) )
form4 = update.formula(form2, . ~ . - plm::lag(Immig_pc, 1) + plm::lag(Imm_mus_pc, 1) )
form5 = update.formula(form2, . ~ . - plm::lag(Immig_pc, 1) + plm::lag(Imm_nonEU_pc, 1) )
form6 = update.formula(form2, . ~ . + plm::lag(Imm_sal_2, 0) )

trust.mods.we = lapply(list(form1, form2, form3, form4, form5, form6), 
                       plm, data=mig.plm.we, model="within", effect="indiv")
sapply(trust.mods.we, function(x) round(pwartest(x, order=1)[["p.value"]][[1]], 3) )
sapply(trust.mods.we, function(x) round(pcdtest(x, test="cd")[["p.value"]][[1]], 3) )
sapply(trust.mods.we, function(x) round(sigma(x), 3) )

vcm = lapply(trust.mods.we, function(x) plm::vcovSCC(x, cluster="time")) # cross-sectional heteroskedasticity and correlation
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(trust.mods.we)) {pval[[i]] = pt(abs(coef(trust.mods.we[[i]]) / rse[[i]]), df=400, lower.tail=FALSE) * 2}
screenreg(trust.mods.we, override.se=rse, override.pvalues=pval, digits=3, single.row=FALSE, leading.zero=FALSE, 
          stars=c(0.05))

vars = c("Institutional trust t-1", "Institutional trust t-2", "% foreign-born t-1", "Immigration rate t-1", 
         "Delta GDP growth per capita t0", "GDP growth per capita t-1", "Unemployment rate t-1", 
         "Corruption t-1", "Liberal democracy t-1", "Far right seat share t-1", "Immigrant integration policy t-1", 
         "Muslim immigration rate t-1", "Non-EU immigration rate", "Concern about immigration t0")
htmlreg(trust.mods.we, file="table_S1.htm", override.se=rse, override.pvalues=pval, digits=3,
        single.row=TRUE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
        caption="Immigration and institutional trust, Western Europe, 1-way FE mods, DKSEs")
texreg(trust.mods.we, file="table_S1.txt", override.se=rse, override.pvalues=pval, digits=3,
       single.row=FALSE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
       caption="Immigration and institutional trust, Western Europe",
       booktabs=TRUE, dcolumn=TRUE,)


## Satis as DV, dynamic FE

form1 = diff(Satis, lag=1) ~ plm::lag(Satis, 1:2) + plm::lag(Forborn_stck_pc, 1) + plm::lag(Immig_pc, 1) 
form2 = diff(Satis, lag=1) ~ plm::lag(Satis, 1:2) + diff(GDP_grwth_WB, lag=1) + plm::lag(GDP_grwth_WB, 1) + 
  plm::lag(Unempl_WB, 1) + plm::lag(Corrup_VD, 1) + plm::lag(Libdem_VD, 1) + plm::lag(Forborn_stck_pc, 1) + 
  plm::lag(Immig_pc, 1)
form3 = update.formula(form2, . ~ . + plm::lag(FRight_seat_pr, 1) + plm::lag(BMI_i_ext, 1) )
form4 = update.formula(form2, . ~ . - plm::lag(Immig_pc, 1) + plm::lag(Imm_mus_pc, 1) )
form5 = update.formula(form2, . ~ . - plm::lag(Immig_pc, 1) + plm::lag(Imm_nonEU_pc, 1) )
form6 = update.formula(form2, . ~ . + plm::lag(Imm_sal_2, 0) )

satis.mods.we = lapply(list(form1, form2, form3, form4, form5, form6), 
                       plm, data=mig.plm.we, model="within", effect="indiv")
sapply(satis.mods.we, function(x) round(pwartest(x, order=1)[["p.value"]][[1]], 3) )
sapply(satis.mods.we, function(x) round(pcdtest(x, test="cd")[["p.value"]][[1]], 3) )
sapply(satis.mods.we, function(x) round(sigma(x), 3) )
vcm = lapply(satis.mods.we, function(x) plm::vcovSCC(x, cluster="time")) # cross-sectional heteroskedasticity and correlation
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(satis.mods.we)) {pval[[i]] = pt(abs(coef(satis.mods.we[[i]]) / rse[[i]]), df=400, lower.tail=FALSE) * 2}
screenreg(satis.mods.we, override.se=rse, override.pvalues=pval, digits=3, single.row=FALSE, leading.zero=FALSE, 
          stars=c(0.05))

vars = c("Satisfaction with democracy t-1", "Satisfaction with democracy t-2", "% foreign-born t-1", "Immigration rate t-1", 
         "Delta GDP growth per capita t0", "GDP growth per capita t-1", "Unemployment rate t-1", 
         "Corruption t-1", "Liberal democracy t-1", "Far right seat share t-1", "Immigrant integration policy t-1", 
         "Muslim immigration rate t-1", "Non-EU immigration rate", "Concern about immigration t0")
htmlreg(satis.mods.we, file="table_S2.htm", override.se=rse, override.pvalues=pval, digits=3, 
        single.row=TRUE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
        caption="Immigration and satisfaction with democracy,Western Europe, 1-way FE mods, DKSEs")
texreg(satis.mods.we, file="table_S2.txt", override.se=rse, override.pvalues=pval, digits=3,
       single.row=FALSE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
       caption="Immigration and satisfaction with democracy, Western Europe",
       booktabs=TRUE, dcolumn=TRUE,)


## Democratic mood as DV, dynamic FE

form1 = diff(SupDem, lag=1) ~ plm::lag(SupDem, 1:2) + plm::lag(Forborn_stck_pc, 1) + plm::lag(Immig_pc, 1) 
form2 = diff(SupDem, lag=1) ~ plm::lag(SupDem, 1:2) + diff(GDP_grwth_WB, lag=1) + plm::lag(GDP_grwth_WB, 1) + 
  plm::lag(Unempl_WB, 1) + plm::lag(Corrup_VD, 1) + plm::lag(Libdem_VD, 1) + plm::lag(Forborn_stck_pc, 1) + 
  plm::lag(Immig_pc, 1)
form3 = update.formula(form2, . ~ . + plm::lag(FRight_seat_pr, 1) + plm::lag(BMI_i_ext, 1) )
form4 = update.formula(form2, . ~ . - plm::lag(Immig_pc, 1) + plm::lag(Imm_mus_pc, 1) )
form5 = update.formula(form2, . ~ . - plm::lag(Immig_pc, 1) + plm::lag(Imm_nonEU_pc, 1) )
form6 = update.formula(form2, . ~ . + plm::lag(Imm_sal_2, 0) )

mood.mods.we = lapply(list(form1, form2, form3, form4, form5, form6), 
                      plm, data=mig.plm.we, model="within", effect="indiv")
sapply(mood.mods.we, function(x) round(pwartest(x, order=1)[["p.value"]][[1]], 3) )
sapply(mood.mods.we, function(x) round(pcdtest(x, test="cd")[["p.value"]][[1]], 3) )
sapply(mood.mods.we, function(x) round(sigma(x), 3) )
vcm = lapply(mood.mods.we, function(x) plm::vcovSCC(x, cluster="time")) # cross-sectional heteroskedasticity and correlation
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(mood.mods.we)) {pval[[i]] = pt(abs(coef(mood.mods.we[[i]]) / rse[[i]]), df=400, lower.tail=FALSE) * 2}
screenreg(mood.mods.we, override.se=rse, override.pvalues=pval, digits=3, single.row=FALSE, leading.zero=FALSE, 
          stars=c(0.05))

vars = c("Democratic support t-1", "Democratic support t-2", "% foreign-born t-1", "Immigration rate t-1", 
         "Delta GDP growth per capita t0", "GDP growth per capita t-1", "Unemployment rate t-1", 
         "Corruption t-1", "Liberal democracy t-1", "Far right seat share t-1", "Immigrant integration policy t-1", 
         "Muslim immigration rate t-1", "Non-EU immigration rate", "Concern about immigration t0")
htmlreg(mood.mods.we, file="table_S3.htm", override.se=rse, override.pvalues=pval, digits=3, 
        single.row=TRUE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
        caption="Immigration and support for democracy, Western Europe, 1-way FE mods, DKSEs")
texreg(mood.mods.we, file="table_S3.txt", override.se=rse, override.pvalues=pval, digits=3,
       single.row=FALSE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
       caption="Immigration and democratic support, Western Europe",
       booktabs=TRUE, dcolumn=TRUE,)


### Simulated Effects Plots 
### Note: These project future immigrant stock based on current immigration flows

## 1. Trust

# immigrant stock forecast model
imm.stck.mod.1 = plm(Forborn_stck_pc ~ plm::lag(Forborn_stck_pc, 1) + plm::lag(Immig_pc, 1), 
                     data=mig.plm, model="within", effect="indiv")

mod = trust.mods.ip[[2]] # includes % foreign born and immigration rate
dem.mod = imm.stck.mod.1
round(coef(mod), 2)
n.coef =  length(coef(mod))
n.burn = 100
n.treat = 30
n.yrs = n.burn + n.treat

summary(mig.trim$Immig_within)
sd(mig.trim$Immig_within, na.rm=TRUE)
summary(mig.trim$Forborn_stck_within)
sd(mig.trim$Forborn_stck_within, na.rm=TRUE)

chg.ip = sd(mig.trim$Immig_within, na.rm=TRUE)
pred.data = data.frame(Trust_m1=0, Trust_m2=0, 
                       diffGDP=0, GDPgr=0, 
                       Unemp=0, Corr=0, LD=0,
                       Forborn_stck_pc=0,
                       IP=c(rep(0, n.burn), rep(chg.ip, n.yrs-n.burn)), 
                       ChgTrust=NA, Trust=NA
                       )
n.sims = 1e3
sim.data = array(as.matrix(pred.data), dim=c(dim(pred.data), n.sims))
mod.sims = mvrnorm(n=n.sims, mu=coef(mod), Sigma=plm::vcovSCC(mod, cluster="time"), empirical=TRUE)
dem.mod.sims = mvrnorm(n=n.sims, mu=coef(dem.mod), Sigma=plm::vcovSCC(dem.mod, cluster="time"), empirical=TRUE)
dem.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)
yhat = matrix(NA, ncol=n.yrs, nrow=n.sims)
trust.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)

for(i in 1:(n.yrs-1)) {
  for(k in 1:n.sims) {
    dem.hat[k, i] = dem.mod.sims[k,] %*% as.matrix(sim.data[i, c(8, 9), k])
    sim.data[i+1, match("Forborn_stck_pc", names(pred.data)), k] = dem.hat[k, i] + rnorm(1, 0, sigma(dem.mod))
    sim.data[i, match("ChgTrust", names(pred.data)), k] = sum(sim.data[i, 1:n.coef, k] * coef(mod))
    sim.data[i, match("Trust", names(pred.data)), k] = sim.data[i, match("Trust_m1", names(pred.data)), k] + 
      sim.data[i, match("ChgTrust", names(pred.data)), k]
    sim.data[i+1, match("Trust_m1", names(pred.data)), k] = sim.data[i, match("Trust", names(pred.data)), k]
    sim.data[i+1, match("Trust_m2", names(pred.data)), k] = sim.data[i, match("Trust_m1", names(pred.data)), k]    
    sim.data[i, match("ChgTrust", names(pred.data)), k] = mod.sims[k,] %*% (as.matrix(sim.data[i, 1:n.coef, k]))
    trust.hat[k, i] = sim.data[i, match("Trust_m1", names(pred.data)), k] + sim.data[i,match("ChgTrust", names(pred.data)), k]
  }
}

pe = apply(trust.hat, 2, mean, na.rm=TRUE)
pe.sd = apply(trust.hat, 2, sd, na.rm=TRUE)
u95 = pe + 1.96 * pe.sd
l95 = pe - 1.96 * pe.sd

yr.st = 101
yr1.off = pe[yr.st]
pe = pe-yr1.off
u95 = u95-yr1.off
l95 = l95-yr1.off

col2rgb(plotcol7[3]) / 255

pdf("Fig_2a.pdf", height=4, width=4)
par(mfrow=c(1,1), mar=c(3, 3, 1, .5), tcl=-0.25, cex=1)
plot(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], type="n", ylim=c(-0.5, 0.5), 
     xaxs="i", yaxs="i", xaxt="n", yaxt="n", xlab="", ylab="")
abline(v=seq(10, 30, by=10), lty=3, col=rgb(0, 0, 0, .4))
abline(h=c(-0.25, 0, 0.25), lty=3, col=rgb(0, 0, 0, .4))
abline(v=0, lty=2, col=rgb(0, 0, 0, .75))
axis(side=1, at=seq(100, 130, by=10), labels=seq(0, 30, by=10), cex.axis=0.9, las=1, mgp=c(0.7, 0.2, 0))
axis(side=2, at=c(-0.5, -0.25, 0, 0.25, 0.5), cex.axis=0.9, las=0, mgp=c(1, 0.4, 0))
lines(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], col=rgb(.18, .44, .54, 1), lwd=2)
polygon(x=c(yr.st:(n.yrs-1), (n.yrs-1):yr.st), col=rgb(.18, .44, .54, .5), border=NA, 
        y=c(u95[yr.st:(n.yrs-1)], l95[(n.yrs-1):yr.st]))
mtext(side=1, line=1.1, "Years", cex=1)
mtext(side=2, line=1.6, "Institutional trust", cex=1.1)
text(x=101, y=0.42, "Immigration increases by 1 SD in year 0", cex=1, pos=4)
dev.off()


## 2. Satisfaction

mod = satis.mods.ip[[2]]
dem.mod = imm.stck.mod.1
round(coef(mod), 2)
n.coef =  length(coef(mod))
n.burn = 100
n.treat = 30
n.yrs = n.burn + n.treat

chg.ip = sd(mig.trim$Immig_within, na.rm=TRUE)
pred.data = data.frame(Satis_m1=0, Satis_m2=0, 
                       diffGDP=0, GDPgr=0, 
                       Unemp=0, Corr=0, LD=0,
                       Forborn_stck_pc=0,
                       IP=c(rep(0, n.burn), rep(chg.ip, n.yrs-n.burn)), 
                       ChgSat=NA, Satis=NA
                       )
n.sims = 1e3
sim.data = array(as.matrix(pred.data), dim=c(dim(pred.data), n.sims))
mod.sims = mvrnorm(n=n.sims, mu=coef(mod), Sigma=plm::vcovSCC(mod, cluster="time"), empirical=TRUE)
dem.mod.sims = mvrnorm(n=n.sims, mu=coef(dem.mod), Sigma=plm::vcovSCC(dem.mod, cluster="time"), empirical=TRUE)
dem.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)
yhat = matrix(NA, ncol=n.yrs, nrow=n.sims)
sat.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)

for(i in 1:(n.yrs-1)) {
  for(k in 1:n.sims) {
    dem.hat[k, i] = dem.mod.sims[k,] %*% as.matrix(sim.data[i, c(8, 9), k])
    sim.data[i+1, match("Forborn_stck_pc", names(pred.data)), k] = dem.hat[k, i] + rnorm(1, 0, sigma(dem.mod))
    sim.data[i, match("ChgSat", names(pred.data)), k] = sum(sim.data[i, 1:n.coef, k] * coef(mod))
    sim.data[i, match("Satis", names(pred.data)), k] = sim.data[i, match("Satis_m1", names(pred.data)), k] + 
      sim.data[i, match("ChgSat", names(pred.data)), k]
    sim.data[i+1, match("Satis_m1", names(pred.data)), k] = sim.data[i, match("Satis", names(pred.data)), k]
    sim.data[i+1, match("Satis_m2", names(pred.data)), k] = sim.data[i, match("Satis_m1", names(pred.data)), k]    
    sim.data[i, match("ChgSat", names(pred.data)), k] = mod.sims[k,] %*% (as.matrix(sim.data[i, 1:n.coef, k]))
    sat.hat[k, i] = sim.data[i, match("Satis_m1", names(pred.data)), k] + sim.data[i,match("ChgSat", names(pred.data)), k]
  }
}

pe = apply(sat.hat, 2, mean, na.rm=TRUE)
pe.sd = apply(sat.hat, 2, sd, na.rm=TRUE)
u95 = pe + 1.96 * pe.sd
l95 = pe - 1.96 * pe.sd

yr.st = 101
yr1.off = pe[yr.st]
pe = pe-yr1.off
u95 = u95-yr1.off
l95 = l95-yr1.off

col2rgb(plotcol7[2]) / 255

pdf("Fig_2b.pdf", height=4, width=4)
par(mfrow=c(1,1), mar=c(3, 3, 1, .5), tcl=-0.25, cex=1)
plot(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], type="n", ylim=c(-0.5, 0.5), 
     xaxs="i", yaxs="i", xaxt="n", yaxt="n", xlab="", ylab="")
abline(h=c(-0.25, 0, 0.25), lty=3, col=rgb(0, 0, 0, .4))
abline(v=0, lty=2, col=rgb(0, 0, 0, .75))
axis(side=1, at=seq(100, 130, by=10), labels=seq(0, 30, by=10), cex.axis=0.9, las=1, mgp=c(0.7, 0.2, 0))
axis(side=2, at=c(-0.5, -0.25, 0, 0.25, 0.5), cex.axis=0.9, las=0, mgp=c(1, 0.4, 0))
lines(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], col=rgb(.24, .30, .54, 1), lwd=2)
polygon(x=c(yr.st:(n.yrs-1), (n.yrs-1):yr.st), col=rgb(.24, .30, .54, .5), border=NA, 
        y=c(u95[yr.st:(n.yrs-1)], l95[(n.yrs-1):yr.st]))
mtext(side=1, line=1.1, "Years", cex=1)
mtext(side=2, line=1.6, "Satisfaction with democracy", cex=1.1)
text(x=101, y=0.42, "Immigration increases by 1 SD in year 0", cex=1, pos=4)
dev.off()


## 3. Support for democracy

mod = mood.mods.ip[[2]]
dem.mod = imm.stck.mod.1
round(coef(mod), 2)
n.coef =  length(coef(mod))
n.burn = 100
n.treat = 30
n.yrs = n.burn + n.treat

chg.ip = sd(mig.trim$Immig_within, na.rm=TRUE)
pred.data = data.frame(SupDem_m1=0, SupDem_m2=0, 
                       diffGDP=0, GDPgr=0, 
                       Unemp=0, Corr=0, LD=0,
                       Forborn_stck_pc=0,
                       IP=c(rep(0, n.burn), rep(chg.ip, n.yrs-n.burn)), 
                       ChgSupDem=NA, SupDem=NA
)
n.sims = 1e3
sim.data = array(as.matrix(pred.data), dim=c(dim(pred.data), n.sims))
mod.sims = mvrnorm(n=n.sims, mu=coef(mod), Sigma=plm::vcovSCC(mod, cluster="time"), empirical=TRUE)
dem.mod.sims = mvrnorm(n=n.sims, mu=coef(dem.mod), Sigma=plm::vcovSCC(dem.mod, cluster="time"), empirical=TRUE)
dem.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)
yhat = matrix(NA, ncol=n.yrs, nrow=n.sims)
sup.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)

for(i in 1:(n.yrs-1)) {
  for(k in 1:n.sims) {
    dem.hat[k, i] = dem.mod.sims[k,] %*% as.matrix(sim.data[i, c(8, 9), k])
    sim.data[i+1, match("Forborn_stck_pc", names(pred.data)), k] = dem.hat[k, i] + rnorm(1, 0, sigma(dem.mod))
    sim.data[i, match("ChgSupDem", names(pred.data)), k] = sum(sim.data[i, 1:n.coef, k] * coef(mod))
    sim.data[i, match("SupDem", names(pred.data)), k] = sim.data[i, match("SupDem_m1", names(pred.data)), k] + 
      sim.data[i, match("ChgSupDem", names(pred.data)), k]
    sim.data[i+1, match("SupDem_m1", names(pred.data)), k] = sim.data[i, match("SupDem", names(pred.data)), k]
    sim.data[i+1, match("SupDem_m2", names(pred.data)), k] = sim.data[i, match("SupDem_m1", names(pred.data)), k]    
    sim.data[i, match("ChgSupDem", names(pred.data)), k] = mod.sims[k,] %*% (as.matrix(sim.data[i, 1:n.coef, k]))
    sup.hat[k, i] = sim.data[i, match("SupDem_m1", names(pred.data)), k] + sim.data[i,match("ChgSupDem", names(pred.data)), k]
  }
}

pe = apply(sup.hat, 2, mean, na.rm=TRUE)
pe.sd = apply(sup.hat, 2, sd, na.rm=TRUE)
u95 = pe + 1.96 * pe.sd
l95 = pe - 1.96 * pe.sd

yr.st = 101
yr1.off = pe[yr.st]
pe = pe-yr1.off
u95 = u95-yr1.off
l95 = l95-yr1.off

col2rgb(plotcol7[4]) / 255

pdf("Fig_2c.pdf", height=4, width=4)
par(mfrow=c(1,1), mar=c(3, 3, 1, .5), tcl=-0.25, cex=1)
plot(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], type="n", ylim=c(-0.5, 0.5), 
     xaxs="i", yaxs="i", xaxt="n", yaxt="n", xlab="", ylab="")
abline(h=c(-0.25, 0, 0.25), lty=3, col=rgb(0, 0, 0, .4))
abline(v=0, lty=2, col=rgb(0, 0, 0, .75))
axis(side=1, at=seq(100, 130, by=10), labels=seq(0, 30, by=10), cex.axis=0.9, las=1, mgp=c(0.7, 0.2, 0))
axis(side=2, at=c(-0.5, -0.25, 0, 0.25, 0.5), cex.axis=0.9, las=0, mgp=c(1, 0.4, 0))
lines(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], col=rgb(.13, .56, .55, 1), lwd=2)
polygon(x=c(yr.st:(n.yrs-1), (n.yrs-1):yr.st), col=rgb(.13, .56, .55, .5), border=NA, 
        y=c(u95[yr.st:(n.yrs-1)], l95[(n.yrs-1):yr.st]))
mtext(side=1, line=1.1, "Years", cex=1)
mtext(side=2, line=1.6, "Support for democracy", cex=1.1)
text(x=101, y=0.42, "Immigration increases by 1 SD in year 0", cex=1, pos=4)
dev.off()

